#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# SConstruct  (Madagascar Script)
#
# Purpose: Obtain zero offset CRS parameters using VFSA.
#
# Site: https://www.geofisicando.com
# 
# Programmer: Rodolfo A. C. Neves (Dirack) 19/09/2019
#
# Email: rodolfo_profissional@hotmail.com
#
# License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.

__author__="Rodolfo Dirack <rodolfo_profissional@hotmail.com>"
__version__="1.0"
__site__="https://www.geofisicando.com"

Help('''
    CRS parameters optimization using VFSA
    ======================================

    Obtain RN, RNIP and BETA parameters using Very Fast Simulated Annealing (VFSA)
    global optimization

    Author: %s
    Version: %s
    Site: %s
    ''' % (__author__,__version__,__site__))

# Madagascar library
from rsf.proj import *

# Python math library
import math

# Ploting functions
def wiggle(title):
    return '''
    wiggle wheretitle=top yreverse=y transp=y label1=Time unit1=s label2=CMP unit2=km pclip=99.5 title="%s" 
    ''' % title

def grey(title):
    return '''
    grey label1=Time unit1=s label2=CMP unit2=km pclip=99.5 title="%s" 
    ''' % title

def grey3(title):
	return '''
	byte |
	transp plane=23 |
	grey3 flat=n frame1=500 frame3=80 frame2=200
	label1=Time unit1=s 
	label3=Offset unit3=km 
	label2=CMP unit2=km
	title="%s" point1=0.8 point2=0.8 
	''' % title

Help('''
### Kirchhoff modeling ###

Use Kirchhoff modeling in a Gaussian reflector linear velocity model
velocity increases with depth with a 0.5 velocity gradient
''')

Flow('gaussianReflector',None,
     '''
     math d1=0.01 n1=2001 o1=-5 unit1=km label1=Offset
     output="4-3*exp(-(x1-5)^2/9)"
     ''')

# Define a gaussian velocity Model
Flow('velocityModel','gaussianReflector',
     '''
     window min1=0 max1=10 |
     spray axis=1 n=451 d=0.01 o=0 label=Depth unit=km |
     math output="1.5+0.5*x1+0.0*x2"
     ''')

Flow('reflectorDip','gaussianReflector','math output="2/3*(x1-5)*input" ')

# Kirchhoff Modeling
Flow('dataCube','gaussianReflector reflectorDip',
     '''
     kirmod cmp=y dip=${SOURCES[1]} 
     nh=161 dh=0.025 h0=0
     ns=401 ds=0.025 s0=0
     freq=10 dt=0.004 nt=1001
     vel=1.5 gradz=0.5 gradx=0.0 verb=y |
     put d2=0.0125 label3="CMP" unit3="Km" label2="Offset" unit2="Km" label1=Time unit1=s
     ''')

m0=5
t0=1.1
v0=1.5

Help('''
### Very Fast Simulated Annealing (VFSA) ###

Obtain RN, RNIP and BETA parameters for m0=%g t0=%g v0=%g
using non-hyperbolic CRS aproximation and VFSA
''' % (m0,t0,v0))

Flow('crsParameters','dataCube',
	'''
	vfsacrsnh om0=%g v0=%g ot0=%g verb=y repeat=10
	''' % (m0,v0,t0))

Help('''
### Calculate approximation error and plot ###

Compare approximation surface with exact traveltime extracted from seismic data
''')

Flow('dataReflectionSurface','dataCube','envelope | max1 | window n1=1 | real')
Plot('dataReflectionSurface',
	'''
	grey color=j bias=2 scalebar=y barlabel=Time barunit=s barreverse=y 
	title="Modeled traveltime surface" label1=Half-Offset unit1=km label2=Midpoint unit2=km
	''')
Flow('crsAppSurface',['dataReflectionSurface','crsParameters'],
	'''
	nhcrssurf param=${SOURCES[1]} m0=%g v0=%g t0=%g verb=y
	''' % (m0,v0,t0))
Plot('crsAppSurface',
	'''
	grey color=j bias=2 scalebar=y barlabel=Time barunit=s barreverse=y 
	title="Non-hyperbolic CRS m0=5Km" label1=Half-Offset unit1=km label2=Midpoint unit2=km
	''')

Flow('error',['crsAppSurface','dataReflectionSurface'],
	'''
	add scale=1,-1 ${SOURCES[1]} | 
	math output="abs(input)" 
	''')
Plot('error',
	'''
	grey color=j scalebar=y barlabel=Time barunit=s barreverse=y 
	title="Approximation error" label1=Half-Offset unit1=km label2=Midpoint unit2=km
	maxval=1 minval=0
	''')

Result('errorAndCRSSurfaces',['dataReflectionSurface','crsAppSurface','error'],'SideBySideAniso')

End()

